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1. Motivations and Objectives . . 

The problem of time discontinuity (or jump condition) in the coalescence/dispersion 
(C/D) mixing model is addressed in first part (PDF) of this work. A C/D mix- 
ing model continuous in time is introduced. With the continuous mixing model, 
the process of chemical reaction can be fuUy coupled with mixing. In the case 
of homogeneous turbulence decay, the new model predicts a pdf very close to a 
Gaussian distribution, with finite higher moments also close to that of a Gaussian 
distribution. Results from the continuous mixing model are compared with both 
experimental data and numerical results from conventional C/D models. 

The effect of Coriolis forces on compressible homogeneous turbulence is studied 
using direct numerical simulation (DNS). The numerical method used in this stu y 
is an eigth order compact difference scheme. Contrary to the conclusions reached 
by previous DNS studies on incompressible isotropic turbulence, the present results 
show that the Coriolis force increases the dissipation rate of turbulent kinetic energy, 
and that anisotropy develops as the Coriolis force increases. The Taylor-Proudman 
theory does apply since the derivatives in the direction of the rotation axis vanishes 
rapidly A closer analysis reveals that the dissipation rate of the incompressible 
component of the turbulent kinetic energy indeed decreases with a higher rotation 
rate, consistent with incompressible flow simulations (Bardina ), while the dissipa- 
tion rate of the Compressible part increases; the net gain is positive. Inertial waves 
axe observed in the simulation results. 

2. Work Accomplished 

2.1 PDF Turbulence Model for Combustion 
In Collaboration with J.Y. Chen 

Accurate prediction of turbulent reacting flows requires the solution of an evo- 
lution equation for the probability density function (pdf) of the thermo-chemical 
variables using Monte Carlo simulation. Since the pdf equation, like most equations 
describing turbulent motion, is not closed, closure models have to be devised. For 
the pdf of scalars, the terms in the pdf equation that need modeling are molecular 
mixing and turbulent convection. The present work deals with the modeling of 

molecular mixing. . ' , . 

Most of the mixing models are based on the coalescence/ dispersion (C/D) model 
by Curl 6 This model is known to have deficiencies, and efforts had been made to 
correct these deficiencies, for example, Janicka, et al. 1979 and Pope 1982. The 
most recent efforts have been devoted to the problem of coupling between mixing 
and chemical reaction. Chen and Kollmann 4 proposed a reaction conditioned model 
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that allows correct prediction of combustion in the flame-sheet regime. Norris and 
Pope 10 proposed a new model based on ordered pairing that aimed at the same end. 

All the existing models suffer in one respect, namely, they are discontinuous in 
time: once a pair of particles are chosen to participate in mixing, their properties will 
jump abruptly regardless of the step size of the time integration. This phenomenon, 
clearly non-physical, could cause difficulty in coupling the processes of mixing and 
fhpmiral reaction. In the present work, a new model that is continuous in time is 
proposed. With this new model, the processes of molecular mixing and chemical 
reaction can be fully coupled. 

In the case of homogeneous turbulence decay of a scalar, one expects a Gaussian 
distribution for the pdf, and finite values for the higher moments. Pope 11 pointed 
out the modified Curl model could not produce the correct pdf for this problem, 
and the higher even moments from that model tend to infinity; he suggested an age. 
biased sampling process to overcome these shortcomings. The present continuous 
model, as we will show, predicts a pdf distribution very close to Gaussian for ho- 
mogeneous turbulence decay, and gives finite higher moments with values close to 
that of a Gaussian distribution. 

The continuous mixing model is applied to the study of both non-reacting and 
reacting flows, and the results are compared with earlier calculations by Hsu 7 as 
well as with experimental data. 

2.1.1 Molecular Mixing Model 

The evolution equation of a single point probability density function of scalar ran- 
dom variables ipi , ..., ip n — representing the species mass fraction and temperature 
— can be written as 

N 

pd t P + pv a d Q P + p^d^ i, ..., iPn)P} 

t=i 

= -d a (p < v'^\(f>k(x) = tp k > P) 

N N 

~/>ZZ 81,*, (< = 'h>P) 

i= 1 3=1 

where the terms represent the rate of time change, mean convection, chemical re- 
action, turbulent convection, and molecular mixing, respectively; P is the density- 
weighted joint pdf: 

P = pP/fi, 

e is the scalar dissipation: 

€ij = Dd a <pid a <pj, 

(where D is the diffusion coefficient), and < x\y > denotes the mathematical ex- 
pectation of a random function x conditioned upon y. 

The left hand side of the above equation can be evaluated exactly and requires 
no modeling; the right hand side terms contain the conditional expectation of the 
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velocity fluctuation and the conditional expectation of the scalar dissipation, which 
arenew unknowns and require modeling. In the present work we concen rate on 
the modeling of the second term, namely, the conditional expectation of 
dissipation, referred to as molecular mixing in the Mowing. 

The Modified Curl Model 

The simplest and most used mixing model is the modified Curl model, wic^ 
assumes binary interaction between sample fluid particles. As describe y P 
in a Monte Carlo simulation, the continuous pdf is replaced by delta functions 


P*{t, t) = jz £ W ~ &*(<)), 


where each delta function represents one sample fluid particle of an ensemble of N 
particles. The evolution of P* entails the movement of the particles m the rp space, 

or the evolution of the individual values of <fi n s. * . . 

With the modified Curl model, the change of <f> n due to molecuiar mmng 
achieved by the following binary interaction process: divide the flow domai 
small ceUs, each containing N sample particles. Given a smaU time interval St and 
a turbulent time scale r, select randomly N mx pairs of particles, 

Nmx = 

{C = 6.0) and let a pair, say, m and n, mix as follows 

4>n(t + fo) ~ A<f> m (t) + (1 — A)<f> n (t) ^ 

<t>m(t + fit) = A<p n (t) + (1 — A)(f>m(t) ^) 

where A = 0.5£, with £ a random variable uniformly distributed on the interval 
[0,1]. The remaining N - 2 N mx particles remain unchanged: 

4> n (t + ^0 = 0n(O 

This model does not represent the true physical process since the properties of 
the sample particles change discontinuously regardless the size of the time mtervid 
St. This deficiency can be best illustrated by rearranging eq. (1) and dividing it by 

(Mt + St)-Mt)) _ - <f,„{t)) (5) 

St ot 

The derivative ^ does not exist because as St goes to zero the right hand side of 
the equation becomes infinite since both A and the difference between ^(t) and 
<t, n (t) are finite. This means that there is a sudden jump in the value of the scalar 
quantities, which is typical of a Poisson process, but is non-physical in the presen 
case since the flow properties of turbulence are continuous. 
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Continuous Mixing Model 

One can see from the previous section that the modified Curl model relies on the 
parameter N mx to control the extent of mixing. On the individual particle level, 
it assumes complete mixing once the particle is selected as one of the mixing pair, 
without considering the size of St. 

In order to achieve continuous mixing, we propose the following model: during a 
time interval St , we assume that all the particles within a cell participate in mixing. 
The extent of the mixing is controlled at the individual particle level. That is to 
say, the N particles within a given cell are randomly grouped into N/2 pairs; the 
properties of all the particles change according to eqs. (1) and (2) The extent 
of mixing now has to be controlled at the individual particle level through the 
parameter A, which is redefined as 


where C‘ = 2.0. With this new definition, eq. (3) can be written, in the limit 
St — 0, 

= c4(<Mi) - «<))• 

at t 

The above equation states that the change of <p n due to mixing is proportional to 
the difference between <p m and 4> n , and inversely proportional to the turbulence 
time scale r. 

The Coupling of Mixing and Reaction 

The processes of mixing and chemical reaction are essentially decoupled when one 
uses the discontinuous C/D models. In contrast, with the above continuous model, 
coupling becomes natural since, for a given particle, mixing and chemical reaction 
can be described with a single equation: 

at r 

where w n is the chemical source term. 

Since the continuous mixing model allows full coupling of the reaction and mixing 
processes, the C/D model with reaction zone conditioning by Chen and Kollmann 4 
can be easily implemented in the present model to simulate the fast reaction in the 
flame sheet regime. Here a modified finite difference version of eq. 12 has to be 
used since w n is infinity in case of fast reaction. 

2.1.2 Results and Discussions 

The continuous mixing model described in the previous section has been vali- 
dated using both non-reacting and reacting flow test cases. The results and their 
comparisons with earlier calculations (Hsu 7 ) using the modified Curl model as well 
as with experimental data are presented in this section. 
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Homogeneous Turbulence Decay of a Scalar 

The case of decaying fluctuation of a passive scalar in homogeneous turbulence 
is used to test the continuous mixing model. The initial condition is 




N / 2 N 

^26{rp- 1 )+ ^2 + 

n=JV/2+ 1 


n=l 


that is, in the Monte Carlo simulation, half of the particles are ascribed the value 

1, and the other half -1. . , . . ■ 

The pdf distribution of the normalized variable (ip- < <P >)/<*> where < <p > 
the mean and a is the standard deviation, in the homogeneous turbulence decay 
problem converges to a single curve after certain time, and the correct distribution 
should be Gaussian. Fig. 1 and 2 are the pdf distributions from the modified Curl 
model and the present model, both compared to a normal distribution. One can 
see that the pdf from the modified Curl model deviates considerable from Gaussian, 
while the result from the present model is fairly close to a Gaussian distribution 
The evolution history of the rms and fourth and sixth moments of the scalar 
fluctuation are calculated using both the modified Curl model and the present 
continuous model. Fig. 3 shows the results from the modified Curl model. One 
can see that although the rms from that model behaves well, the fourth and sixth 
moments grow quickly out of bound, oscillating at a level several order of magnitudes 
higher than the value of Gaussian distribution. These results are similar to what 
Pope 11 had observed. Fig. 4 shows the results for the same set of quantities from the 
present model. The rms behaves similar to that from the modified Curl model. The 
fourth and sixth moments, on the other hand, are quite different from those of the 
previous model; they rise smoothly to the value predicted by Gaussian distribution. 
Although the values do not seem to converge, they remain finite, and are of the 
same order of magnitude as that of the Gaussian distribution. 

The above results clearly demonstrated the advantage of the present model over 
that of the modified Curl model. Pope 11 had devised an age biased scheme that 
achieved the same end, which required an additional variable, namely the age ot 
the particles, and two extra adjustable parameters. In contrast, the present mode 
needs no extra work or parameters. 


Heated Turbulent Jet . , , . , 

Extensive experimental results for a heated turbulent plane jet have been reported 

by many authors (Bashir, et al. 2 , Browne et al. 3 , Uberoi and Singh , Jenkins ^ An- 
tonia, et al. 1 ). The turbulent jet has a slightly higher temperature than the ambient. 
Measurements of both the the mean temperature and the rms of temperature fluc- 
tuations were given. We compared the solutions for the heated turbulent jet from 
the new model with experimental data as well as with previous solutions (Hsu ) 

obtained using the modified Curl model. 

In the present study, a combined CFD-Monte Carlo algorithm is used. The 
mean flow field is obtained by solving the Navier-Stokes equation and a two-equation 
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turbulence model using a finite difference scheme. The temperature is treated as a 
conserved scalar and simulated by the pdf equation. 

Fig. 5 shows the comparisons of the the mean temperature distribution from 
the pdf Monte Carlo simulations and experimental data from various authors. The 
figure shows that both mixing models predict the mean temperature distribution 
accurately. 

The results for standard variation, or rms, of the temperature distribution are 
given in Fig. 6. Although the two solutions do not show significant difference, the 
new model seems to agree slightly better with the experimental data. The skewness 
and the flatness, i.e., the third and fourth moments of the temperature fluctuation, 
are given in Figs. 7 and 8; the comparisons show that in the present case of a 
turbulent jet the statistical behavior of the new model is similar to that of the 
modified Curl model. 

Hydrogen-Fluorine Diffusion Flame 

The continuous model has been applied to the H 2 — F 2 diffusion flame studied 
earlier by Hsu 7 . The flow conditions are set according to an experiment performed 
by Hermanson and Dimotakis 6 . The flame consists two streams. The upper stream 
contains §6% of N 2 and 4% of F 2 , the flow velocity is U\ = 22 m/s ; the lower stream 
contains 96% of N 2 and 4% of H 2 , with velocity U 2 = 8.8 m/ s. The estimated 
Damkohler number ranged from 25 to 130 (Hermanson and Dimotakis 6 ), and a fast 
chemistry model is deemed appropriate in the calculation. Again a modified version 
of eq. 12 is used to accommodate the fast chemistry. 

Fig. 9 shows the temperature rise due to combustion. In the figure, St is the 
shear layer thickness determined by 1% of the temperature rise, AT is the actual 
temperature rise due to combustion, ( the two streams have the same temperatures 
initially,) and A T ad f is the adiabatic flame temperature assuming complete reaction. 
Details on the flow conditions can be found in Hermanson and Dimotakis 6 . The 
agreement between numerical predictions and experimental data is fairly good, and 
a comparison of the results from the continuous model and that from the modified 
Curl model shows that both performed well for this case. 

Combination with Reaction Zone Conditioning 

Chen and Kollmann 4 developed a mixing model based on reaction zone condi- 
tioning, aimed at the coupling of reaction and mixing. We have shown in Section 2.3 
that with the present model, the processes of reaction and mixing can be fully cou- 
pled; therefore it is only natural to apply the reaction zone conditioning suggested 
by Chen and Kollmann here. 

The H 2 — F 2 diffusion flame problem is reformulated such that the chemical reac- 
tion is confined to a very narrow zone near stoichiometry. By applying reaction zone 
conditioning to the continuous mixing model, we were able to produce a scatter plot 
of the temperature vs. mixture fraction in which all the points reached the equi- 
librium temperature. This result is shown in Fig. 10. The mixture fraction here is 
defined as the molar concentration of fuel divided by the total molar concentration, 
and stoichiometry is located at / = 0.5. 
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2.1.3 Conclusions 

A turbulence mixing model that is continuous in time has been introduced. The 
deficiency of non-physical jump condition in the mixing process is removed m the 
new model. It has been shown that the new model is supenor to the existing 
modified Curl model (Janicka, et al. 8 ) in that it can predict a Gaussian distribution 
and finite higher moments in the case of homogeneous turbulence decay; it has 
accomplished what the age biased sampling scheme (Pope ) is designed for without 
the extra parameters required by that scheme. The numerical results from the 
present model compare well with experimental data. 


2.2 DNS of Homogeneous Compressible Turbulence in a Rotating frame 

Rotation is an important factor in many flow phenomena in nature and in engi- 
neering. Problems that are strongly affected by rotation include flows in turboma- 
chinery, large scale motions in the atmosphere and oceans, and galactic motions. 
Turbulence is important in all the above examples. In order to model turbulence 
in a rotating frame, we need a better understanding on the effects of rotation on 

Several experiments had been carried out by various researchers 14 . However, 
the conclusions are inconsistent. Some results show that rotation increases the 
turbulence dissipation rate, others suggest the opposite. 

Numerical studies of turbulence in a rotating frame are few. In fact the only 
two are incompressible flow studies, by Bardina, et al . 17 and by Speziale et . , 
using the same computer program. Studies other than the present one dealing with 
compressible flows in a rotating frame are not known to the authors. 

Bardina et al . 17 performed both large eddy and full simulations of incompress- 
ible isotropic turbulence. Their results show that, in the case of incompressible 
turbulence, the dissipation rate decreases with increasing rotation rate, and that 
anisotropy does not develop as a result of the Coriolis force. Speziale et al. con- 
firmed these results using a smaller Rossby number, i.e. a faster rotation. 

In general, compressible turbulence in a rotating frame is not homogenous or 
isotropic because of the existence of centrifugal force and density fluctuations. How- 
ever when the rotation rate is low, or when the solution domain is very close to the 
rotation axis, the centrifugal force can be neglected and the flow is approximately 
homogeneous. In the present study, our goal is to identify the effects of the Cori- 
olis force, as opposed to the effects of the centrifugal force or the combined effects 
of the centrifugal and Coriolis forces. Therefore, the centrifugal force is dropped 
from the governing equations, regardless of whether or not it is negligible. A future 
study will concentrate on the non-homogenous effects of the centrifugal force and 

the combined effects of the two forces. . .. 

The numerical results presented in the following show that although the Coriolis 
force does not appear explicitly in the turbulence kinetic energy equation, its effect 
can be significant. Anisotropy does develop as a result of the Coriolis force. The 
dissipation rate of the compressible and incompressible modes of the turbulent ki- 
netic energy behave differently under Coriolis forces, and the combined effect is a 
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increased dissipation rate. 

2.2.1 Method 

The governing equations for compressible flows in a rotating frame are 
P,t + (puj),j = 0 

(pUi) t t + = P,i + T ij,j ~ P^ij k^klm ft j ft — 2/?e t jfcftjUfc 

(pe),t + (puje) tj = -pujj + u i jT ij - (kTj) d 
p = pRT 

The last two terms in the second equation represent, respectively, the effects of 
centrifugal force and the Coriolis force. The centrifugal force is a nonuniform body 
force which would induce nonuniform pressure and density distributions; its effect 
on turbulence is to destroy homogeneity. The Coriolis force, on the other hand, is 
a uniform body force, provided that the velocity field is uniform. In the case of 
incompressible flows, the centrifugal force can be included into a modified pressure 
term and thus does not appear explicitly; however, this can not be done in the case 
of compressible flows. 

In order to identify the homogeneous effect of the Coriolis force alone, and to 
compare results with incompressible flow simulations, we simply take away the 
centrifugal force in our study. This, of course, greatly simplifies the problem: With 
a homogeneous initial condition, the flow field remains homogeneous under the 
Coriolis force. The simplified problem is a hypothetical one and can not normally 
be found in nature except in very restricted situations such as the ones mentioned 
in the introduction. 

The governing equations are solved using a compact difference scheme. Lele 19 has 
shown that such schemes have spectral accuracy and are well suited for DNS. The 
spatial derivatives are approximated by 8th order finite differences, and the time 
integration is performed using a three-stage Runge-Kutta time marching scheme. 

The initial condition for the rotating flow simulation is obtained from a simulation 
of isotropic decaying turbulence. The initial condition for the decaying turbulence 
is generated using a given 3D spectrum of the form 



After the decaying turbulence shows the correct skewness, the Coriolis force is 
imposed onto the flow field. 

2.2.2 Results and Discussion 

As a primary study, three cases have been calculated, with rotation in the z- 
direction: The nondimensional rotation rates are ft = 0, 40, 80, which correspond 
to turbulent Rossby numbers of oo, 0.0075, and 0.00375. These Rossby numbers 
are much smaller than those considered in References 17 and 18. 
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Figures 11 and 12 show the development of a ID two-poirit correlation and a ID 
spectrum in the turbulent flow field, where square symbols are initial conditions an 
fines represent later times. Figure 13 shows a comparison of the energy spectrum 
obtained using two different grids, one double the size of the other. The agreemen 


affirms sufficient resolution. 

Figures 14-19 show the time history of various statistics of the three cases, where 
time is nondimensionafized using the initial eddy turn over time. 

Figure 14 is the ratio of the turbulent kinetic energy, where ql is the initial 
turbulent kinetic energy and q* the value at later times. It is immediately clear 
that the Coriolis force causes a faster decay of the turbulence, a conclusion that is 
contrary to what was observed in the incompressible flow simulations by previous 
researchers 17 ’ 18 . Nevertheless, as we shall see in the following, the present results 
are in fact consistent with those of Refs 4 and 5. 

To prove that the present results are consistent with the previous incompressible 
simulations, we look at the ratio of the means of the divergence squared and vorticity 
squared, which are proportional to, respectively, the compressible and incompress- 
ible turbulent kinetic energy. Figures 15 and 16 show that while the dissipation 
rate of the compressible part of the turbulent kinetic energy increases with larger 
Coriolis forces, the dissipation rate of the incompressible part decreases. The com- 
bined effect, as we have seen from Figure 14, is an increase in the total dissipation 
rate. Since the simulations of References 4 and 5 are for incompressible flows, it is 
not surprising that only a decrease in dissipation rate was observed. . 

In both References 17 and 18, no anisotropy was observed. Bardina et al ' 
speculated that a reorganization to two-dimensional flows could occur at ig er 
rotation rates, while Speziale et al. suggested that Taylor-Proudman reorganization 
would not occur in a rapidly rotating isotropic turbulence. The present results show 
definite signs of anisotropy as soon as the Coriolis force is turned on. The anisotropy 
can be clearly observed from the time history of the normal components of the 
Reynolds stress, given in Figures 17, 18, and 19 for the three rotation rates. With 
the Coriolis force, the x- and y-components of the normal Reynolds stress decrease 
much faster than the z-component. In fact, it seems that the increased dissipation 
is primarily in the x- and y-directions, and the z-component of the turbulent kinetic 
energy is not affected much. 

The Taylor-Proudman theory stipulates that under strong rotation, with the 
rotation axis parallel to, say, the z axis, then for any Sow variable u , we have 
(du/dz) = 0. Figures 20 and 21 are the contours of the x-component of velocity 
plotted on three planes. Without the Coriolis force (Figure 20), the flow struc- 
ture remains isotropic. With the Coriolis force, columns soon appear in the flow 
field suggesting that a Taylor-Proudman reorganization has occurred. Using larger 
Rossby numbers such as those used in References 17 and 18, the results show no 
appreciable anisotropy, hence there appears to be no contradiction. 

Finally, inertial waves can be observed from the time history of various statistics 
presented in Figures 14-19. The frequency is approximately fl/ir, the theoretical 
intrinsic frequency for rotating flows. 


25 


A. T. Hsu 


3. Future Work 

3.1 PDF 

The pdf model is being implemented into an existing compressible flow solver: 

the RPLUS2D code. Work on implementing pdf to RPLUS3D will follow. 

3.2 DNS 

To understand the effect of compressibility, incompressible homogeneous turbu- 
lence in a rotating frame will be studied. 
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Figure 2. Asymptotic pdf distribution for a scalar in homo- 
geneous turbulence. — present model; - - - Gaureian. 



Figure 3. Evolution of moments from tbs modifie d Curl 

t ivmM — standard deviation, 0-01 x fourth central 

moment, - - 0.0001 x «ixth central moment, o 0-01 x fourth 
nw r 1 ^ * for Gtutciin distribution, O 0.0001 x sxxth moment 
for Gioinu distribution. 



Figure 4. Evolution of moments from the press* model — 
standard deviation, • • • 0*1 x fourth central moment, - — 
0X1 x sixth central moment, o 0.1 x fourth m o men t for 
Gaussian distribution, □ 0.01 x sixth moment for Gaussian 
distribution. 
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5 Mean temperature in heated plane jet. — contin- 
uoua mixing model, - - - modified Curi modd, A Browne ei 
al. ( 19 M) 7 ° Baahir and Uberoi (1975), o Uberoi and Smgh 
(1975), • Jenldna and Goldadimidt (1973). 




Figure 6 . RMS of temperature variance in heated plane jet. 
_ eontinuoM mitring modd, - - - modified Curl modd, o 
Antonia et ai (1983), o Baahir and Uheroi (1975), A Uberoi 
and Singh (1975). 



Figure 7. Sh twne w of temperature ▼ariance in heated pUae 
continuous ^tiring model, • • - modified Curi nwid. 




Fijnire 9 . Temperature rise in »n Fj 
I^ZLoum i£dng model, * - - modified Curl model, o 



r.dfc ^ 

Mixtur* Fraction. / 


Figure 10. Joint pdf between temperature and mixture frac- 
tionon the centerline of the H r F t diffnrion lUme; reacUon 
U restricted to a narrow rone at rtoichiometry; 
obtained by applying reaction rone condiUonmg (Cben and 
Kollmann, 1990) to the continuour mixing modd. 
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